Proteomic insight into soybean response to flooding stress reveals changes in energy metabolism and cell wall modifications

Soybean is a legume crop enriched with proteins and oil. It is frequently exposed to anthropogenic and natural flooding that limits its growth and yield. Current study applied gel-free proteomic techniques to unravel soybean response mechanism to flooding stress. Two-days-old soybeans were flooded for 4 days continuously and root samples were collected at days 2 to 6 for proteomic and enzymatic analyses. Age-matched untreated soybeans were collected as control. After protein extraction, purification and tryptic digestion, the peptides were analyzed on nano-liquid chromatography-mass spectrometry. A total of 539 and 472 proteins with matched peptides 2 or more were identified in control and flooded seedlings, respectively. Among these 364 proteins were commonly identified in both control and flooded soybeans. Fourty-two protein’s abundances were changed 4-fold after 2-days of flooding stress as compared to starting point. The cluster analysis showed that highly increased proteins included cupin family proteins, enolase, pectin methylesterase inhibitor, glyoxalase II, alcohol dehydrogenase and aldolase. The enzyme assay of enolase and pectin methylesterase inhibitor confirmed protein abundance changes. These findings suggest that soybean adopts the less energy consuming strategies and brings biochemical and structural changes in the cell wall to effectively respond to flooding stress and for the survival.


Introduction
Soybean (Glycine max (L.) Merr.) is an important legume that is enriched with proteins and oil contents [1]. Frequent flooding due to climatic changes and ill-drained fields is one of the abiotic stresses that reduce its growth and yield [2]. Flooding initially causes damage to the roots [3], reduces the nutrient uptake [4] and decreases the nitrogen fixation capacity [5]. Flooding stress reduces biomass, tap-root length, and pod number, inhibits carbon/nitrogen content in root/nodule, decreases nodule dry weight, and grain yield in soybean [6]. These reports suggest that flooding is a major constraint on growth and yield of soybean. Root is an important primary organ to feel the effects of flooding stress. Flooding reduces the root dry weight first [7]. Oxygen transport from the air to the roots is important for root physiology [8]. Flooding causes oxygen deficiency leading to hypoxia or anoxia as oxygen moves ten thousand times slower in water than in the air [8,9]. Plants respond to flooding stress by formation of adventitious roots [7,10] and aerenchyma formation [7]. Adventitious roots formation benefits the plant growth during flooding exposure [11]. Flooding stress did not affect root growth of submergence-tolerant rice genotypes [12]. Roots undergo structural and functional alterations at the cellular, molecular and phenotypic level to deal with the flooding stress [13]. Roots rapidly use starch reserves for limiting the damage and maintaining the growth [3].
Proteomic techniques found extensive applications in investigating effects of flooding stress and flooding stress-responsive proteins. Proteins belonging to the categories of glycolysis, fermentation, detoxification of reactive oxygen species, anaerobic catabolism, storage, stress, development, cell organization, transport, signaling and amino acid metabolism-related proteins were changed in abundance under flooding stress [14][15][16]. Proteins related to the cell wall lignification were suppressed [17]. Protein abundances of energy-related proteins were raised whereas those involved in protein folding and cell structure organization were lowered in flooded soybean [18]. Kamal et al. [19] reported a decrease in sucrose metabolism-related proteins but increase in fermentation-related proteins in soybean cotyledon under flooding stress. Photosynthesis, RNA, DNA, signaling, and the tricarboxylic acid cycle were changed in abundance leaf, hypocotyl and root of soybean under flooding stress [20]. Proteomics approaches have also been applied on subcellular level to reveal localized cellular responses and investigate communications among subcellular components during flooding stress. In the plasma membrane, proteins related to signaling, stress and the antioxidative system were increased; whereas, reactive-oxygen species scavenging enzymes activities were retarded in the cell wall [21]. Protein metabolism-related proteins were decreased in the nucleus and also proteins related to electron transport chain were suppressed in the mitochondria [21]. The soybean responses to flooding stress are being studied at various levels utilizing proteomic approaches. Current proteomic study was designed to analyze response mechanism of soybean to continuous four days flooding stress.

Protein extraction
An amount of 500 mg of root was weighed on balance and ground under liquid nitrogen using a mortar and pestle. The powder was transferred to pre-cooled (4˚C) acetone solution containing 10% trichloroacetic acid and 0.07% 2-mercaptoethanol. The mixture was vortexed and sonicated for 10 min. The suspension was incubated for 1 h at -20˚C and then centrifuged at 9,000×g at 4˚C for 20 min. The pellet was washed twice with 0.07% 2-mercaptoethanol in pre-cooled acetone and dried. It was resuspended in lysis buffer (7 M urea, 2 M thiourea, 5% CHAPS, 2 mM tributylphosphine) by vortexing for 1 h at 25˚C and centrifuged at 25˚C with 20,000×g for 20 min. The supernatant was collected as protein extract. Bovine serum albumin was used as standard for protein concentration calculations through Bradford assay (0 to 2 mg/mL standard curve range) [22].

Protein purification and digestion for mass spectrometry analysis
Protein extracts of 100 μg were purified with methanol and chloroform to remove detergent from the samples. For purification and digestion of extracted proteins, methodology described by Khan and Komatsu [23] was followed. Briefly, 150 μL protein extract was mixed with 600 μL methanol. The resulting suspension was mixed with 150 μL chloroform and 450 μL water and centrifuged at 20,000×g for 10 min to achieve phase separation. The upper aqueous phase was discarded and 450 μL methanol was added slowly to the lower phase. The samples were further centrifuged at 20,000×g for 10 min, and the obtained pellets were dried at room temperature. The dried samples were reduced with 50 mM dithiothreitol for 30 min at 56˚C, followed by alkylation with 50 mM iodoacetamide for 30 min at 37˚C in the dark. Alkylated proteins were digested with trypsin and lysyl endopeptidase at 1:100 enzyme/protein concentrations at 37˚C for 16 h. The resulting tryptic peptides were acidified in 20% formate in milli-Q water and analyzed by nano-liquid chromatography (LC) mass spectrometry (MS).

Nanoliquid chromatography-tandem mass spectrometry analysis
A nanospray LTQ Orbitrap mass spectrometer (Thermo Fisher Scientific, San Jose, CA, USA) was operated in data-dependent acquisition mode with the installed XCalibur software (version 2.0.7, Thermo Fisher Scientific). The nanoLC-MS conditions and method as described by Khan and Komatsu [23] was followed. Peptides in 0.1% formic acid were loaded onto a C18 PepMap trap column (300 μm ID × 5 mm; Dionex, Germering, Germany) of an Ultimate 3000 NanoLC system. The peptides were eluted from the trap column with a linear acetonitrile gradient (8-30% over 120 min) in 0.1% formic acid at a flow rate of 200 nL/min. The eluted peptides were separated and sprayed onto a C18 capillary tip column (75 μm ID × 120 mm; Nikkyo Technos, Tokyo, Japan) at a spray voltage of 1.5 kV. Full-scan mass spectra were acquired in the nanospray LTQ Orbitrap MS system over 400-1500 m/z with a resolution of 30,000. A lock mass function was used for high mass accuracy [24]. The six most intense precursor ions were selected for collision-induced fragmentation in the linear ion trap at a normalized collision energy of 35%. Dynamic exclusion was employed within 90 s to prevent the repetitive selection of peptides [25].

Protein identification by Mascot search
Proteins were identified from a soybean peptide database (73,320 sequences, 29,844,971 amino acid residues) constructed from the soybean genome database (Phytozome version 9.1, http:// www.phytozome.net/soybean) [26] using the Mascot search engine (version 2.4.0.2, Matrix Science, London, UK) through Proteome Discoverer software (version 2.3.2, Thermo Fisher Scientific). For the Mascot searches, the carbamidomethylation of cysteine was set as a fixed modification and the oxidation of methionine was set as a variable modification. Trypsin was specified as the proteolytic enzyme and one missed cleavage was allowed. Peptide mass tolerance was set at 10 ppm, fragment mass tolerance was set at 0.8 Da, and peptide charge was set at +2, +3, and +4. An automatic decoy database search was performed as part of the search. Mascot results were filtered with the Mascot Percolator package to improve the accuracy and sensitivity of peptide identification. False discovery rates for peptide identification of all searches were less than 1.0%. False discovery rates for peptide identification of all searches were less than 1%. Peptides with > 13 (p < 0.05) percolator ion score were used for protein identification.

Differential analysis of acquired mass spectrometry data
The Mascot results were exported for SIEVE software analysis (version 2.1; Thermo Fisher Scientific). SIEVE compares the relative abundances of peptides and proteins between control and experimental groups. For analysis, chromatographic peaks detected by MS were aligned and peptide peaks were detected as a frame on all parent ions scanned by MS/MS using 5 min of frame time width and 10 ppm of frame m/z width. Area of chromatographic peak within a frame was compared for each sample and ratios between samples were determined for each frame. The frames with MS/MS scan were matched to the Mascot search results. The ratio of peptides between samples was determined from the variance-weighted average of the ratios in frames, which matched the peptides in the MS/MS spectrum. The ratios of peptides were further integrated to determine the ratios of corresponding proteins. Total ion current was used for normalization of differential analysis of protein abundance. The outliers of ratio were deleted in frame table filter based on the frame area. The minimum requirement for protein identification was two matched peptides. Significance of protein abundance between samples was analyzed (p < 0.05).

Cluster and in silico protein-protein interaction analyses
Protein ratios obtained from SIEVE software analysis were subjected to cluster analysis using Genesis software (version. 1.8.1; http://genome.tugraz.at) [27]. Genesis software was downloaded from the mentioned website and license was obtained on request. Cluster analysis was performed using hierarchical clustering with a Euclidean distance metric and a centroid linkage clustering method. The clustered proteins alignment in treatment was used for heat map generation in control. Clustered proteins were analyzed for in silico protein-protein interactions utilizing online STRING database (version 11.0; https://string-db.org).

Functional categorization
The functional categories of identified proteins were determined through MapMan bin codes using MapMan software (http://mapman.gabipd.org) [28]. Log2-fold change values with Glyma ID codes were used for the MapMan analysis. Proteins were categorized based on their BIN names and codes as per MapMan functional categories (1-35).

Analysis of enzyme activities
Enolase. A quantity of 200 mg of root was homogenized in lysis buffer (20 mM Tris-HCl pH 7.5, 1 mM EDTA, 1 mM 2-mercaptoethanol) prepared in milli-Q water. The suspension was centrifugation at 20,000×g at 4˚C for 30 min. The supernatant was collected as enzyme extract. Protein concentrations were estimated by Bradford assay using bovine serum albumin as standard (0 to 2 mg/mL standard curve range) [22]. A reaction mixture consisting of 100 mM triethanolamine (pH 7.4), 120 mM KCl, 2.25 mM 2-phosphoglycerate, 0.2 mM 2-NADH, 30 mM MgSO 4 , 1.75 mM ADP, 10 units pyruvate kinase, and 15 units L-lactic dehydrogenase was used for enzymatic assay. Enzyme extract of 100 μL was mixed with 900 μL of reaction mixture and vortexed. The absorbance was measured in 1 cm standard cuvette at 340 nm using a UV/Vis spectrophotometer (DU Series 700, Beckman Coulter, CA, USA) [29,30].
Plant invertase/pectin methylesterase inhibitor superfamily. Plant invertase assay was performed by slightly modifying protocol of Huang et al. [31]. The extraction procedure was performed on ice. A weight of 200 mg of soybean roots was used for enzyme extraction. Roots were ground into fine powder in liquid nitrogen and extracted in buffer that consisted of 50 mM HEPES-KOH, pH 7.4, containing 5% Polyvinyl pyrrolidone, 1 mM EDTA, 1 mM EGTA, 1 mM PMSF, 5 mM DTT, 0.1% Triton X-100, and 1% glycerol. The homogenate was centrifuged at 15000×g for 20 min at 4˚C. The supernatant was collected as the enzyme crude extract. The crude extract was vacuum-filtered through bottle-top vacuum filters (pore size: 0.45 μm). The filtrate was used for enzyme assay. An enzyme extract of 100 μL was mixed with 900 μL of reaction mixture and reduction in absorbance was measured at 340 nm in a 1 cm standard cuvette using a UV/Vis spectrophotometer (DU Series 700, Beckman Coulter, CA, USA).

Statistical analysis
Enolase and Pectin methylesterase activities were analyzed for statistical significance through One-way ANOVA. The Post Hoc Duncan's multiple range test was applied for analyzing significances at specific sample points through SPSS (version 21.0). The statistical differences were represented with different alphabets.

Identified proteins in soybean root under flooding stress
To identify differentially changed proteins in soybean root, a gel-free proteomic technique was used to analyze the protein profiles of soybeans that had been flooded continuously for 4 days. A total of 539 and 472 proteins with matched peptides 2 or more were identified in control (S1 Table) and flooding-stressed soybean roots (S2 Table), respectively. Out of the total identified proteins, 364 were commonly identified in control and flooding-stressed plants (S3 Table; Fig  2). Among these 364 proteins, protein abundances of 42 proteins were changed 4-fold in flooding-stressed plants after 2-days of flooding ( Table 1). The mass spectrometry proteomics data files were deposited to the ProteomeXchange Consortium via the iProX partner repository [32].

Identified proteins belonged to diverse functional categories
The total identified proteins in control (539) and flooded soybean (472) had 364 commonly changed proteins. The total identified proteins were functionally categorized according to MapMan codes (Fig 3A). Maximum number belonged to protein-metabolism category with 152 in control and 117 in flooded soybeans. The second major category was stress-related proteins with 33 identified in control and 34 in flooded seedlings. The other differentially changed
In control plant proteins aligned against flooded cluster II, abundances of RmlC like cupins superfamily protein, formate dehydrogenase, and urease accessory protein G were very less as compared to age-matched flooded plants. In control plant proteins aligned against flooded cluster III, LEA domain containing protein, cupin family protein, stress induced protein, embryonic cell protein 63, RmlC like cupins superfamily protein, annexin 2 and beta amylase 5 were decreased in abundance throughout the growth period; whereas, these proteins were increased in flooded plants.

Enolase and plant invertase/pectin methylesterase inhibitor show highly significant response to flooding stress
The enzyme enolase (Glyma16g32960) which is also called phosphopyruvate hydratase is an important enzyme of glycolysis, was analyzed for activity changes under flooding stress. The protein abundance of enolase was highly increased under initial 2 days of flooding stress (10.28) and decreased gradually latter at day 3 and 4 of flooding stress (5.82 & 2.84) (Fig 6A). While in control plants, there was no appreciable elevation with increasing age. The results of enolase activity assay followed the pattern of protein abundance changes. The enzyme activity was highly raised from first to second day of flooding (160.65 to 720.15 unit/mg protein) and gradually d-lowered at days 3 and 4 of flooding (600.25 & 470.58 unit/mg protein, respectively) (Fig 6B). The changes in activity were significant as compared to those observed in control plants and also among the different flooding durations.
Plant invertase also called pectin methylesterase inhibitor (PMEI) showed a high increase in protein abundance (Fig 7A). The protein abundance raised from 1.57 after 1 day of flooding towards maximum of 49.65 at the end of 3 days flooding. It was reduced at the end of 4 day of flooding to a level of 28.87. The enzyme activity of PMEI was analyzed in control and flooded plants (Fig 7B). PMEI activity gradually elevated from 90.77 after 1 day flooding to a highest of 390.47 unit/mg protein at the end of 4-days flooding period. The activity changes were statistically significant in the last 2 days of flooding.

Discussion
Flooding stress causes injury in the soybean [16]. In the current study, continuous flooding stress was applied to the soybeans for 4 days and protein abundance changes were analyzed through gel-free proteomic technique. The study was conducted to unravel the mechanism involved in soybean responses to continuous flooding stress. Flooding stress brought huge abundance changes in many physiologically important proteins. Among the functionally

PLOS ONE
important proteins, abundances of cupin family protein, RmlC like cupins superfamily protein, enolase, plant invertase/pectin methylesterase inhibitor protein, Arabidopsis thaliana peroxygenase 2, seed maturation protein, glyoxalase II 3, alcohol dehydrogenase 1 and aldolase supefamily protein were significantly increased under flooding stress as compared to starting point 2(0) as well as control plants (Table 1). Proteins related to protein metabolism categories such as synthesis, degradation and folding were also raised in abundance under flooding stress from 4.06 to 10-times as compared to starting point (Table 1). Protein metabolism-related proteins such as ribosomal proteins belonging to different protein families were among the highest interacting proteins when analyzed by STRING.
RmlC-like cupin superfamily proteins and cupin family proteins, which include storage proteins belonging to the development category, were highly increased in abundance under flooding stress. Cupin are functionally very diverse family of proteins [33] and play role in seedling development in soybean [34]. Cupins and seed maturation proteins with nutrient reservoir activity, are development-related storage proteins that were also previously reported to be increased in flooded soybean roots possibly due to delayed degradation [35,36]. The results of the current study suggest delayed use of cupins as storage proteins in the initial 3 days of flooding stress as against control plants where their abundance was quite low. The other types of cupins modify the structure of cell wall as phosphomannose isomerase modifies mannose derivatives [37]. Cupins such as dTDP-rhamnose enzymes produce activated rhamnose as germin cross-link the plant cell-wall components [38,39]. Hence cupins are vital for cell survival through modification of cell wall. The increased abundance of cupins in the flooded soybean may point out towards their role in maintaining cell wall integrity under flooding stress.
Glyoxalase II was increased in flooded 7-fold as compared to starting point and 3-fold as compared to 4-days age-matched control. This enzyme is involved in detoxification of methylglyoxal whose production is increased many-folds under abiotic stress [40]. Methylglyoxal II is produced as by-product of metabolic pathways such as glycolysis and from photosynthesis intermediates (glyceraldehyde-3-phosphate & dihydroxyacetone phosphate). Methylglyoxal is a reactive cytotoxin that can cause lipid peroxidation, oxidation of proteins & fatty acids and disruption of membranes [41,42]. Methylglyoxal is detoxified by glyoxalase system consisting of glyoxalase I and glyoxalase II that catalyze conversion of methylglyoxal to D-lactate while using glutathione as co-factor [40]. The increased protein abundance of glyoxalase II in current study showed an increase in detoxification of methylglyoxal as a defense effort by soybean.
Aldolase superfamily protein abundance was increased at 2 nd and 3 rd days of flooding as compared to control plants. Aldolase is an enzyme that brings conversion of fructose bisphosphate to glyceraldehyde-3-phosphate and dihydroxyacetone phospate, an important step of glycolysis. The enzyme is also involved in gluconeogenesis and calvin cycle [43,44]. Nuclear isoform of fructose-bisphosphate aldolase regulates expression of its own gene as well as other genes by acting as DNA-binding protein [45]. Aldolase is induced under hypoxia that may result from abiotic stress [46]. Aldolase is linked with tonoplast for the activity of V-ATPase in salt-stressed Mesembryanthemum crystallinum that results in sodium ion accumulation in vacuole as a defense strategy [47]. Fructose bisphosphate aldolase is speculated in integration of signals linked to the growth, development, and sugar anabolism [48]. In soybean exposed to flooding stress, aldolase protein abundance was increased [49]. Fructose bisphosphate aldolase is induced by various abiotic stresses in Arabidopsis [50]. The enzyme is also involved in plant development, metabolism and abiotic stress responses [51]. In the current study, increased protein abundance of aldolase depicts increased rate of glycolysis under flooding stress as plant had limited means to generate energy due to blockage of oxidative phosphorylation.
In the current study, protein abundance of the enolase was increased under flooding stress. The enzyme activity changes also followed the pattern of increase. The enzyme enolase which is also called phosphopyruvate hydratase is an important enzyme of glycolysis, responsible for conversion of 2-phosphoglycerate to phosphoenol pyruvate that ultimately leads to pyruvate formation along-with energy generation. This particular enolase is localized in the cytosol, binds magnesium ion and possesses phosphopyruvate hydratase activity. Moreover, it is involved in vacuole fusion [26]. Enolase is induced in maize under anaerobic conditions [52]. Enolase has also been shown linked to the tonoplast for enabling V-ATPase activity [47]. Increase in enolase abundance has been reported in soybean facing flooding stress [49,53]. The results of the present study are in agreement with previous reports indicating that enolase as glycolytic enzyme might have helped in increasing frequency of glycolysis for generating energy under flooding stress.
Alcohol dehydrogenase 1 protein abundance was highly increased under flooding stress as compared to age-matched control plants. Under anaerobic conditions such as flooding, plants ferment glucose to ethanol in the presence of alcohol dehydrogenase. Fermentation thus produces small amount of ATP for life continuity along-with glycolysis [54]. Proteomic and transcript abundances of alcohol dehydrogenase are highly increased in soybean under flooding stress [36,49,55]. Activities of alcohol dehydrogenase were remarkably increased in soybean leaf under flooding stress [20]. From the previous reports as well as results of current study, the evidence of alcohol dehydrogenase induction and shifting of metabolism to anaerobic mode is confirmed. Soybean used anaerobic fermentation to increase its ATP for survival under flooding stress.
Plant invertase/pectin methylesterase inhibitor was increased in protein abundance and activity. The enzyme activity was much higher when measured at the end of 3 rd and 4 th day of flooding stress. Pectin plays roles in controlling cell wall porosity [56], cell adhesion [57] and a key factor in plant development [58,59]. Pectin methylesterase (PME) brings esterification. The extent of methylesterification determines the susceptibility of the plant cell wall to the pectin-degrading enzymes [60]. Plant PME activity generates methanol as a signal of the damaged self, leading to regulate the transcription of pathogen-related PME inhibitor (PMEI) genes [61]. Studies suggest that inhibitory activities of PMEIs are crucial depending on the cell wall environment and different specificities for target PMEs for ensuring a development-and/or stress-dependent adjustments in cell wall [62]. Plant invertase/PMEI abundance and/or activity increased in soybean under flooding stress in current study as well as previous findings by Oh and Komatsu [49] and Yasmeen et al. [53]. These reports suggest that cell wall brings readjustments in its structure and mechanics as a mechanism to deal with the flooding stress.

Conclusions
Flooding acts as abiotic stress for soybean that brings hypoxic or anoxic conditions on the plant. Soybeans respond to flooding stress by altering its basic metabolic modes. It restricts the normal metabolism and brings reduction in ATP yielding and high energy consuming processes. Plant accelerates glycolysis as glycolytic enzymes such as aldolase, enolase etc. increase their protein abundances and activities. Side-wise, after glycolysis, pyruvate undergoes fermentation pathway to yield ethyl alcohol. Multi-faceted Cupins and toxics scavenging glyoxalases also play crucial roles in flooding stress responses. Cell wall being outer boundary of plant cell is at high exposure to flooding stress but brings alterations and rearrangements in its structure and mechanics through various enzymes such as pectin methylesterase inhibitors to cope with the flooding stress. Thus, soybean brings biochemical and structural changes to effectively respond to flooding stress and adopts the less energy consuming strategies for the survival.
Supporting information S1